Quantum dynamics in splitting a harmonically trapped Bose-Einstein condensate 
by an optical lattice: Truncated Wigner approximation 
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We study the splitting of a harmonically trapped atomic Bose-Einstein condensate when we con- 
tinuously turn up an optical lattice (or a double-well) potential. As the lattice height is increased, 
quantum fluctuations of atoms are enhanced. The resulting nonequilibrium dynamics of the frag- 
mentation process of the condensate, the loss of the phase coherence of atoms along the lattice, and 
the reduced atom number fluctuations in individual lattice sites are stochastically studied within 
the truncated Wigner approximation. We perform a detailed study of the effects of temperature 
and lattice height on atom dynamics, and investigate the validity of the classical Gross-Pitaevskii 
equation in optical lattices. We find the atom number squeezing to saturate in deep lattices due 
to nonadiabaticity in turning up of the lattice potential that is challenging to avoid in experiments 
when the occupation number of the lattice sites is large, making it difficult to produce strongly 
number squeezed (or the Mott insulator) states with large filling factors. We also investigate some 
general numerical properties of the truncated Wigner approximation. 

PACS numbers: 03.75.Lm,03.75.Kk,03.75.Gg 



I. INTRODUCTION 

The experimental progress in loading ultra-cold atomic 
gases in weakly coupled mesoscopic traps formed by peri- 
odic optical potentials has provided a dilute atomic sys- 
tem with strongly enhanced interactions. Examples of 
this progress include experiments on the Bose-Einstein 
condensate (BEC) coherence [1-5], the atom number- 
squeezed states [6], the Mott insulator (MI) phase tran- 
sition [7-11], atom dynamics [11-15], and fermionic sys- 
tems [16-18]. Ultra-cold atoms trapped in an optical lat- 
tice resemble traditional condensed matter crystal lattice 
systems [19], but optical lattices are amenable to much 
higher experimental control. The lattice strength can be 
easily modified [2, 6-9, 20] and it could be possible to en- 
gineer, e.g., complex lattice geometries [21], interatomic 
interactions [22, 23], and the spatial profiles of the atomic 
hopping amplitude between the lattice sites [24-29]. 

In order to produce atomic lattice systems close to 
their ground state, atoms need to be adiabatically loaded 
into the optical lattices by continuously turning up the 
lattice potential. The BEC is initially confined in a har- 
monic trap. The increasing strength of the lattice poten- 
tial reduces the tunneling amplitude between the neigh- 
boring lattice sites and the system becomes more strongly 
interacting. The strong interactions enhance quantum 
fluctuations, eventually destroying the long-range phase 
coherence of the atoms and fragmenting the BEC. If the 
turning-up of the lattice potential is not adiabatic, the 
quantum fluctuations of the atoms can be far from their 
ground state properties, resulting in complex dynamics. 

In order to preserve adiabaticity during the turning- 
up of the optical lattice potential, the rate of change of 
the Hamiltonian has to be slower than any time scale of 
the system. The required slow ramping-up time gen- 
erally makes it very difficult to reach the MI ground 



state experimentally, when there are many atoms per 
lattice site, as we demonstrated in Ref. [30]. The sit- 
uation is different from the MI state experiments [7-11] 
with small filling factors, because of rapidly emerging en- 
ergy gap and larger quantum fluctuations in the case of 
small atom numbers. The nonadiabaticity of the turning 
up of the optical lattice has recently been theoretically 
studied also by using classical Gross-Pitaevskii equation 
(GPE) (without thermal and quantum fluctuations) : The 
buildup of an inhomogeneous phase profile [31, 32] was 
identified in the classical analysis as a potential signature 
for nonadiabaticity. 

In this paper we study matter wave dynamics beyond 
the classical mean-field theory by considering a harmoni- 
cally trapped BEC that is dynamically split by an optical 
lattice potential. In our numerical model the quantum 
and thermal fluctuations of the atoms are included in a 
classical stochastic field description within the truncated 
Wigner approximation (TWA). We provide a more de- 
tailed and extended analysis of our previous results [30] 
by also studying the limits of validity of the classical GPE 
and the general numerical properties of the TWA. The 
lattice systems we consider have large occupation num- 
bers per site. Recent experimental studies by Orzel et 
al. [6] on the squeezing of atom number fluctuations have 
been performed in the similar limit of strong optical lat- 
tices and large filling factors. The states with squeezed 
atom number fluctuations with high occupation numbers 
are of great interest in precision measurements, as they 
may allow Heisenberg limited interferometry with poten- 
tial technological applications [33]. 

We study the nonequilibrium dynamics of the BEC 
that is split by the lattice potential, describing the loss of 
phase coherence along the lattice, the atom number fluc- 
tuations within individual lattice sites, and the fragmen- 
tation process of the BEC. We show in numerical TWA 
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simulations how the atom number fluctuations evolve in 
time after the turning up of the lattice potential from 
almost Poisson initial value to a strongly reduced result 
that depends on the parameters of the system. We also 
show that the squeezing of atom number fluctuations sat- 
urates for deep lattices, as a consequence of nonadiabatic 
turning up of the lattice potential. The saturation of the 
number squeezing for deep lattices was experimentally 
observed in Ref. [6]. This system was not tightly elon- 
gated as our ID numerical model, but we can still make 
qualitative comparisons to the experimental data. Al- 
though the saturation was assumed in Ref. [6] to be an 
artifact of the analysis method of the interference mea- 
surement, we also numerically find a similar saturation 
effect in a qualitative agreement with the experiment. 
Since the atom number fluctuations in Ref. [6] were indi- 
rectly detected from the phase noise of the interference 
measurement, not all the experimentally detected phase 
fluctuations may have directly corresponded to the atom 
number squeezing, due to the possibility of nonadiabatic 
loading of the atoms. We find, however, in the TWA sim- 
ulations that considerable atom number squeezing can be 
present even in the nonadiabatic regime in the presence 
of increased phase noise and a nonuniform phase profile. 

In the numerical TWA simulations we first consider a 
BEC that is split by a simple double-well potential and 
compare the numerical results to the phase collapse rate 
calculated for an equilibrium BEC in a double- well poten- 
tial. We then apply the TWA to the splitting of a BEC 
by a periodic optical lattice potential that constitutes the 
main results of the paper. We evaluate the dynamics of 
the phase coherence between atoms occupying different 
lattice sites and the atom number fluctuations in individ- 
ual sites during and after the the splitting of the BEC. 
We study in detail different cases when we vary the final 
height of the lattice potential, the initial temperature and 
nonlincarity. In order to investigate the limits of valid- 
ity of the classical GPE, we also systematically compare 
the TWA simulation results to the GPE results by vary- 
ing the initial temperature and the final lattice height. 
In shallow lattices and at low temperatures the two ap- 
proaches yield very similar results for the coherence prop- 
erties and atom statistics. However, as the lattice height 
and/or the temperature are increased the validity of the 
GPE becomes poorer. 

We also study the generation of the initial state in 
the TWA simulations using an ideal gas. The stochastic 
noise is sampled for noninteracting atoms, before turn- 
ing up the nonlinearity in order to produce the initial 
state of an interacting system. We find that this tech- 
nique works better at low temperatures in which case the 
results are close to the ones obtained by more accurate 
TWA calculations. 

Since the TWA simulations return symmetrically order 
expectations values, we also investigate the importance 
of the particular stochastic representation of the density 
matrix. We demonstrate by numerical examples that the 
correct transformation between symmetric and normal 



operator orderings is very crucial in obtaining the correct 
physical results, especially at low temperatures. 

The TWA was introduced in nonlinear optics to study 
quantum fluctuations [34] and it has provided success- 
ful descriptions for nonlinear optical squeezing [35]. The 
TWA and other closely related classical field methods 
[30, 36-46] have previously also been applied to the 
studies of atomic BEC dynamics. In our earlier stud- 
ies [30, 43] we found that the TWA is particularly use- 
ful to study bosonic atom dynamics in optical lattices 
in the limit where the full multi-mode dynamics, beyond 
the tight-binding approximation, becomes important and 
when the atom filling factor in the lattice is large. In 
Ref. [43] the TWA was able to describe the experimen- 
tally observed dissipative atom dynamics in tightly con- 
fined, shallow ID optical lattices [11] even in a strongly 
fluctuating quantum system. 

In Sec. II we introduce the TWA and the methods 
to extract the normally ordered physical observables in 
the numerical simulations. We also describe the gen- 
eral numerical approach used in the TWA simulations in 
Sec. II C. The splitting of a harmonically trapped BEC 
by an optical potential is studied in Sec. III. We first 
consider a double-well potential in Sec. Ill A and an op- 
tical lattice with a small number of sites in Sec. IIIB. 
The main numerical results of the paper are presented 
in Sec. Ill C where we address an optical lattice with 
a larger number of sites. The effects of the final lat- 
tice height on the atom number fluctuations and matter 
wave coherence are studied in Sec. Ill CI. We analyze 
the adiabaticity of the turning up of the lattice poten- 
tial and present qualitative comparisons to the experi- 
ment of Ref. [6]. The effects of initial temperature and 
nonlinearity on the dynamics are studied in Sec. IIIC2. 
We compare the TWA simulation results to the classical 
GPE results in Sec. Ill C 3 by varying the initial temper- 
ature and the final lattice height. The change of tem- 
perature during the splitting is analyzed in Sec. Ill C 4. 
An alternative method for generating the initial state of 
the TWA simulations by turning up the nonlineariaty is 
studied in Sec. IIIC5. In Sec. IIIC6 we consider the 
effect of the initial state noise sampling in the TWA and 
in Sec. Ill C 7 the importance of symmetric operator or- 
dering in the Wigner representation. A few concluding 
remarks are made in Sec. IV. The Bogoliubov approx- 
imation in a discrete tight-binding approximation is in- 
troduced in Appendix A and the three-body losses are 
estimated in Appendix B. 



II. OVERVIEW OF THE METHOD 

A. Truncated Wigner approximation 

We assume that the bosonic atoms are trapped in a 
tight elongated (along the x direction) cigar-shaped (pro- 
late) trap V 3 d(x, y, z) = m[Lu 2 x 2 + uj 2 ^ 2 + z 2 )]/2, with 
the trap frequencies satisfying to = lo x <C wj_ = u> y = u> z . 
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Here the radial frequency is denoted by lo± and the axial 
frequency by uj. We ignore the density fluctuations along 
the transverse directions in order to obtain an effective 
fD Hamiltonian: 

J dxft(x)ft(x)i>{x)1>(x). (1) 



+ 



Here 



C = T + V h (x) + V (x,t)-n, 



(2) 



includes the kinetic energy f = —h 2 d 2 /(2m), the har- 
monic trapping potential along the axial direction Vh = 
mu) 2 x 2 /2 7 and the time-dependent periodic optical lattice 
potential, or in the case of the double- well a Gaussian po- 
tential barrier, V {x, t). The atom mass and the chemical 
potential are denoted by m and /i, respectively. The ef- 
fective ID interaction strength is given by g\u = 2ftw±a, 
where a is the scattering length. 

We study the dynamics of a finitc-temperaturc BEC 
within the TWA when we load the atoms into an opti- 
cal lattice by ramping up the periodic lattice potential. 
The TWA may be obtained by using the familiar tech- 
niques of quantum optics [47, 48] to derive a generalized 
Fokker-Planck equation for the Wigner distribution of 
the trapped multi-mode BEC [36]. We obtain from the 
Hamiltonian (1): 



di 



/oo ( r 

Jxil — [(c + g 1D {m 2 -l))^ 



1 



■%l>\W(W) + c.c. (3) 



4 S 2 tp5ip : 

Here the density operator of the quantum system is repre- 
sented by a classical quasidistribution function W(ip, tp*) 
of the complex functions that correspond to the 

field operators (-0,-0^). The expectation values in the 
Wigner representation (■ ■ -)w are obtained from the qua- 
sidistribution function: 

(4>*{xi) ■ ■ ■ ip*(x k )ip(x k+ i) ■ ■ ■ ip(xi)) w = 

J d 2 i> w(ip, rw{xi) ■ ■ ■ r(xkWx k+1 ) ■ ■ ■ v(^) . 

(4) 

The expectation values obtained according to the Wigner 
distribution correspond to the expectation values of 
quantum operators that are in symmetric, or Weyl, or- 
der. 

The diffusion matrix of the Fokker-Planck equation 
(3) for W(ip,ip*) vanishes identically and the dynami- 
cal quantum noise acts via third-order derivatives. It is 
useful to write Eq. (3) in terms of stochastic differen- 
tial equations for (ip,ip*) whose ensemble average of the 
dynamics generates the expectation values (4) obtained 



from the quasidistribution function. The TWA consists 
of neglecting the third-order derivatives in Eq. (3), re- 
sulting in a deterministic equation for the classical field 
Tpw which coincides with the Gross-Pitaevskii equation 
(GPE) [36]: 



ih 



dip w (x,t) 
dt 



£ + 9lD(\?Pw{x,t)\ 2 - 1) 1pw(x,t). 



(5) 

Here we have introduced a subscript in ipyy in order to 
emphasize that it denotes the classical Wigner represen- 
tation of the field operator. We have also explicitly in- 
cluded the constant phase factors from Eq. (3). Although 
the time evolution represented by Eq. (5) is deterministic, 
the thermal and quantum fluctuations are still included 
in the initial state of tpw that represents an ensemble of 
Wigner distributed wave functions. The neglected terms 
in the TWA are small when the amplitudes of the Wigner 
distribution are large. It should be emphasized that with- 
out ignoring the third-order derivatives in the Fokker- 
Planck equation no simple stochastic description exists. 
The stochastic representation in the TWA for the field 
operator exists for the states for which W (ipw , tpw) 1S 
positive. The stochastic description is especially useful 
since a single field ipw incorporates both the condensate 
and the noncondensate populations. 

The atoms are initially assumed to be in thermal equi- 
librium in a harmonic trap, before the periodic optical 
lattice potential is turned up. In order to sample the ini- 
tial state stochastically, we solve the quasiparticle excita- 
tions of the BEC within the Bogoliubov approximation. 
We expand the field operator r i)j(x,t = 0) in terms of the 
BEC ground state amplitude do^o, with (djd ) = N , 
and the excited states: 

ip(x) = %1)q(x)o. + ^2 [ u ji x )^j " v j( x )a-l] ■ (6) 

j>0 

Here otj are the corresponding quasiparticle annihilation 
operators, with 



(&%} = nj ee 



1 



exp {ej/ksT) — 1 



(7) 



and ipo is ground state solution of the GPE with the 
chemical potential /x. The quasiparticle mode functions 
Uj (x) and Vj (x) (j > 0) and the corresponding eigenener- 
gies Cj are obtained as solutions to the Bogoliubov equa- 
tions in the subspace that is orthogonal to the ground 
state wave function ip^\ 

(t + 2N Q g 1D \ip a \ 2 'j uj - N a gi D ip 2 Vj = CjUj, 
(t + 2N g 1D \ip 1 2 ) vj - N g 1D ipQ 2 Uj = -ejVj . (8) 

In the TWA we unravel the dynamics into individ- 
ual stochastic trajectories where the classical represen- 
tation ipwix) of the initial state of the quantum field 
operator ip(x, t = 0) in Eq. (6) is sampled according to 
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its Wigner distribution W(ipw, ipw)' so ^ na ^ ^ ne ensem - 
ble average of these individual realizations synthesizes 
the correct quantum statistical correlation functions ac- 
cording to Eq. (4). In particular, in order to construct 
tpw{x,t = 0) from ip(x,t = 0) in Eq. (6), we replace the 
excited state quantum operators (ay, at) (for j > 0) by 
the complex random variables (ctj, a*), obtained by sam- 
pling the corresponding Wigner distribution of the quasi- 
particles. Within the Bogoliubov approximation the op- 
erators (dj,dj) behave as a collection of ideal harmonic 
oscillators whose Wigner distribution in a thermal bath 
may be easily evaluated [48]: 

2 

W(aj,a*) = - tanh (£,) exp [-2|a,| 2 tanh (£,)] , (9) 

where £j = ej/2ksT. The Wigner function is Gaussian 
distributed with the width fij + \ . The nonvanishing con- 
tribution to the width at T = for each mode represents 
the quantum noise. Since the Wigner function returns 
symmetrically ordered expectation values, we have 

{a*a ) w = Jd 2 aj W{ aj , a*)\ aj \ 2 =nj + ^, (10) 

and similarly (aj)w = {a*)w = (ot 2 )w = 0, etc. In 
many realistic physical situations the number of modes 
required to generate the fluctuations in the initial state 
can significantly exceed the lowest energy band in optical 
lattices [43], emphasizing the multi-band nature of the 
TWA. 

We consider large BECs with N » 1. Since the BEC 
is initially weakly-interacting and confined in a harmonic 
trap with no optical lattice, the main contribution to the 
matter wave coherence is due to the thermal and quan- 
tum fluctuations of low-energy phonons that are more im- 
portant than the quantum fluctuations of the initial state 
of the BEC mode. In most cases studied in the paper, we 
sample the quantum noise of the BEC mode according 
to the Wigner distribution of the coherent state [48] : 

W c (a , a*) = | exp (-2|a - iV 1/2 | 2 ) , (11) 

for which (a }w = N^ 2 and {a^a^jw — Nq + |. 
Since we compare the matter wave coherence between 
the atoms in different lattice sites, the global BEC phase 
is unimportant. The advantage of using the coherent 
state description is that the corresponding Wigner dis- 
tribution is positive. As shown later in the paper, the 
assumption of the coherent state for the BEC mode is 
not very important. In fact, we could even treat the 
BEC mode classically without significantly affecting the 
results. In Section III C 6 we compare the sampling of the 
BEC mode according to Eq. (11) to the case in which the 
BEC mode is treated classically having a fixed number 
of atoms. A classical treatment of the ground state does 
not affect the prediction for the phase coherence along 
the lattice, but produces slightly smaller atom number 
fluctuations in the lattice sites. 



Our TWA model for the trapped atomic vapor is for 
the Hamiltonian describing a closed system. In the lat- 
tice experiments the atoms are also coupled to environ- 
ment, resulting in dissipation with the system relaxing 
towards its ground state. We could introduce a more so- 
phisticated model, e.g., by incorporating the three-body 
losses and the spontaneous emission due to the lattice 
lasers. This would introduce also a dynamical noise term 
in Eq. (5). Although the spontaneous emission generates 
an important loss mechanism for the atoms in several 
experimental situations, it can be reduced by further de- 
tuning the lattice laser light from the resonance of the 
atomic transitions. For instance, with intense CO2 lasers, 
that are far off-resonant, the spontaneous emission rate 
can be very low [49]. In Appendix B we estimate the 
importance of the three-body losses of atoms in a ID op- 
tical lattice for a typical set of parameters considered in 
the TWA simulations. 



B. Symmetric ordering 

The Wigner distribution returns symmetrically or- 
dered expectation values for the field operators. This 
means that in the TWA simulations the expectation val- 
ues involving the full multi-mode Wigner fields are sym- 
metrically ordered with respect to every mode. In gen- 
eral, this can significantly complicate the analysis of the 
numerical results from the TWA simulations. For in- 
stance, for the initial state (6) the simple spatial correla- 
tion function in the Wigner representation would actually 
return: 

(r w (x)i>w(x'))w = (ftwfa')) + ^r (x)Mx') 

+ JEk(^k(^')+^(xK(x')], (12) 
1 i=i 

where (■ ■ -}w denotes the expectation value obtained 
from the TWA simulations and (• • •) the normally or- 
dered expectation value of the quantum operators. In or- 
der to extract normally ordered expectation values for the 
correlation functions of the BEC from the TWA simula- 
tions, Steel et al. in Ref. [36] defined a 'condensate mode' 
operator associated with the projection of the stochastic 
field onto the ground state solution. This was used to cal- 
culate the phase diffusion of a spatially static, harmoni- 
cally trapped BEC. Since here we study the splitting of 
a BEC by a periodic optical lattice potential, it is useful 
to define analogously the ground state operators aj for 
each individual lattice site j: 

a j(t) = / dxipQ(x,t)ijtw(x,t) , (13) 

where ipw{x,t) is the stochastic field, determined by 
Eq. (5), and ip (x,t) is the ground state wave function 
at time t, obtained by integrating the CPE in imaginary 
time in the potential V(x,t). The integration is over 
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one lattice site. In the following we also use the same 
description in a double-well potential in which case we 
integrate over the left and the right wells. The impor- 
tance the correct operator ordering is demonstrated by a 
numerical example in Section III C 7. 

With the projection method we can avoid the com- 
plications arising from the symmetrically ordered multi- 
mode field tpw Using the definition in Eq. (13), the nor- 
mally ordered expectation values can be easily obtained 
for each lattice site ground state mode cij. For instance, 
the atom number in the ground state in the jth site reads: 



rij = (a]aj) = {a*a 3 ) w - ^ , 
and the atom number fluctuations in the jth site: 

1/2 



(14) 



Arij 



{{a*ajf) w - (a*jaj)w 



1/2 



(15) 



Similarly, in order to characterize the phase coherence 
along the lattice, we may introduce the normalized first- 
order correlation function Cj between the atoms in the 
central well of the lattice and in its jth neighbor: 



(16) 



In order to obtain directly the fluctuations in the rel- 
ative phase operator ipoj = <f>j — <£>o between the atoms 
in the central lattice site and in its jth nearest neighbor, 
we evaluate Dj defined as: 



D 3 = \ (e^° 



-<Pj)] 



)l- 



(17) 



Here aj = \aj\e Vj and Dj is calculated by normalizing 
a\aj in each stochastic trajectory before the averaging. 
Typically Dj and Cj yield almost equal results and we 
usually only show Cj. 



the occupation number of the central lattice sites is al- 
ways high. For the typical nonlincarity N$gir> = lOOfiuZ, 
with I = (h/muj) 1 / 2 , the corresponding initial Thomas- 
Fermi radius R/l = (3N a g 1D /2hul) 1 / 3 ~ 5.3. In ID the 
strength of interactions is commonly expressed in terms 
of the ratio 7 between the interaction energy and the ki- 
netic energy needed to localize the atoms within the mean 
interatomic distance l/rim, where n\r> is the ID atom 
density [50]. We obtain 7 = mgi£)/(?i 2 ni£>) < 10~ 3 and 
the initial harmonically trapped BEC is well described 
by the GPE and the Bogoliubov theory. 

The time evolution of the ensemble of Wigncr dis- 
tributed wavefunctions, determined by Eq. (5), is unrav- 
eled into stochastic trajectories, where the initial state 
of each realization for the classical stochastic field tpw 
is generated using the expansion (6) with the operators 
replaced by the complex, Gaussian-distributed [Eqs. (9) 
and (11)] random variables (aj,a*j). During the time evo- 
lution we continuously increase the strength of the peri- 
odic optical lattice potential or, in the case of the double- 
well, a repulsive Gaussian potential at the center of the 
trap. The integration of the time dynamics [Eq. (5)] is 
performed using the nonlinear split-step method [51] on a 
spatial grid of up to 4096 points and in several cases with 
the optical lattice potential the sufficient convergence is 
obtained after 600 realizations, whereas in a double-well 
potential typically around 800 realizations are required. 
The convergence is generally slower at higher temper- 
atures. Finite temperature systems require more iter- 
ations and are therefore numerically more demanding. 
This is clear in the sampling of the Wigner distribution 
in Eq. (9) whose width increases with the initial temper- 
ature 7j, indicating higher excited level population and 
more noise in the initial state Wigner distribution. Un- 
like the 3D TWA [38], the ID simulations do not similarly 
depend on the total number of quasiparticle modes and 
we found the calculated results to be unchanged when we 
varied the number of modes. 



III. SPLITTING A BEC BY AN OPTICAL 
POTENTIAL 



C. Numerical approach 



A. Double-well potential 



We study the nonequilibrium quantum dynamics of 
the bosonic atoms within the TWA. The BEC, which 
is initially confined in a harmonic trap, is split by a pe- 
riodic multi-well optical lattice or a double-well poten- 
tial. The atoms are initially assumed to be in thermal 
equilibrium and we first solve the ground state ampli- 
tude profile tpo(x) by evolving the GPE in imaginary time 
in the harmonic trap. We then diagonalize the Bogoli- 
ubov equations [Eq. (8)] in the subspace orthogonal to 
ipo(x) in order to obtain the quasiparticle eigenfunctions 
Uj(x),Vj(x) and the corresponding eigenenergies €j. 

Throughout the paper we use large atom numbers, 
Ao = 2000 atoms at the beginning of the ramp, so that 



Before turning to the dynamical studies of bosonic 
atoms fragmented by an optical lattice, it is useful first 
to introduce the TWA method in a much simpler double- 
well potential case. Since the projection method we 
found very useful in an optical lattice is much less ac- 
curate in a double-well potential for short nonadiabatic 
ramping-up times, we only present here some example 
cases and will publish a more detailed study elsewhere. 

The splitting of a BEC by a double-well potential has 
attracted considerable theoretical interest in the past 
(see, e.g., Refs. [52, 53], and references therein). Also the 
finite-temperature dynamics of atomic BECs in a double- 
well potential has been investigated [54, 55]. Recent ex- 
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periments include, e.g., the use of a double- well BEC as 
a noise thermometer [56]. 

We use the similar set-up as in the early BEC exper- 
iment of Ref. [57], where a harmonically trapped BEC 
was split by a blue-detuned far-off resonant laser beam. 
We model the laser by a Gaussian potential in Eq. (2) 



V (x, t) = A(t) exp(-a; 2 /er) , 



(18) 




cot 



10 



1.0 
0.5 



5 cot 10 



with a — 0.5Z 2 . During the time evolution, we increase 
exponentially the laser potential to some final value Af 
at time r according to A(t) = exp(ret) — 1. In Fig. 1 we 
show the dynamics of the phase coherence C\ between 
the two wells, as defined in Eq. (16). We vary the non- 
linearity and the initial temperature. The ramping-up 
time of the Gaussian potential is lot = 5. As expected, 
both the increase in the nonlincarity and in the initial 
temperature result in a faster reduction of the relative 
phase coherence. However, the projection to the ground 
state in Eq. (13) in the double-well case resulted in no- 
table atom number oscillations and the calculated phase 
coherence values should therefore be considered only as 
an order of magnitude estimates. We can compare the 
zero-temperature case of the collapse time of the relative 
phase coherence between the atoms in the two wells to 
the estimates obtained from the Bogoliubov theory for 
a ground state BEC [58] . If we generalize the results of 
Ref. [58] for the present situation and for arbitrary An, 
we obtain for the collapse time r c : 



T C = 2 4 v3 ln2 (!^r 2/3 41 

V hioi I LoAn 



(19) 



After a time r c , the value of the coherence function C\ 
is reduced to 0.5. The values of n and An in Eq. (19) 
could be obtained from the numerical results of the TWA 
simulations, e.g., by averaging An after the ramping pro- 
cess. If we assume the atom number fluctuations between 
the two wells to be Poissonian (or binomial), we obtain 
for the collapse times lot c = 3.3,2.5,2.1, for the nonlin- 
earities Nogm / (tvjjl) = 100,150,200, respectively. The 
collapse time of the TWA simulations is estimated by de- 
termining the time when C\ = 0.5 and subtracting the 
ramping-up time lot = 5 from this value. For the same 
nonlinearities we obtain lot c > 5, lot c ~ 5, and lot c ~ 4. 

The collapse time for the phase coherence is quite sen- 
sitive to the atom number fluctuations and could be ex- 
perimentally used to measure An. The collapse time in 
the TWA simulations indicates sub-Poisson atom num- 
ber fluctuations. The simple collapse time estimate from 
Eq. (19), that does not take into account the effects 
of the dynamical splitting of the BEC, yields An ~ 
Anpoisson/2. The numerical values of An in the TWA 
simulations are not reliable because of the simple projec- 
tion method used in the double-well problem. 



FIG. 1: The dynamics of the phase coherence during the 
splitting of a BEC by a Gaussian potential barrier. The 
phase coherence between the atoms in the two wells for dif- 
ferent nonlinearities (on left) NogiD / (fuvl) = 100, 150, 200 
(the curves from top to bottom) at the initial temperature 
Ti = 0. The phase coherence at two different initial temper- 
atures (on right) ksTi/huj = and 37.8 (the lower curve) for 
Nogm = 200fkol. In the both cases the final value of the 
Gaussian potential is Af — WOOhui, and its width a = 0.5Z 2 . 
The ramping-up time is lot — 5 and iVo = 2000. 



B. Optical lattice with a small number of sites 

The condensate fragmentation is generally much more 
difficult to analyze when the number of wells is increased 
from two [59]. The TWA, incorporating quantum and 
thermal fluctuations for the complete field operator, be- 
comes a very useful tool in investigating the dynamics 
of the fragmentation process. Unlike in the double-well 
case, the projection method to the lowest mode in each 
site produces stable and more accurate results. 

In this Section we demonstrate the loss of phase coher- 
ence between the atoms in different lattice sites, when the 
BEC is split by continuously turning up an optical lat- 
tice with a small number of sites. We evaluate the phase 
collapse both at T = and at finite temperature. In 
the following Section we turn to a more detailed study of 
the turning up process of the lattice potential when the 
number of lattice sites is larger. 

We start with the system in thermal equilibrium in 
the harmonic trap and continuously turn up an optical 
lattice potential, so that V a (x, t) in Eq. (2) is defined by 



V (x,t) = s(t)E r sm 2 (nx/d) , 
where the lattice photon recoil energy 

E 

Sli V 



2md 2 



(20) 



(21) 



and d = A/2 sin((9/2) denotes the lattice period, obtained 
by two laser beams intersecting at an angle 9. The ad- 
vantage of ID lattices is that the lattice spacing can be 
easily modified by changing the angle between the lasers. 
In the simulations the height of the lattice potential s(t) 
is turned up exponentially during a time r to some final 
value s according to 



s(t) 



1 



(22) 
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FIG. 2: The normalized ground state atom density in a har- 
monic trap without optical lattice (left) and with the lattice 
for s = 37.5. The lattice has clearly pushed the atoms to- 
wards the outer regions of the trap by increasing the radius 
of the atom vapor. The nonlinearity Nogm = 12Qtkol and the 
wavenumber kl — 4. 



for t < t. In this Section we use the lattice spacing d = 
ttI /4 and the nonlinearity N^gm — 120fruil. For instance, 
the ground state calculated at the lattice height s = 37.5 
then occupies ~ 19 sites; see Fig. 2. However, in a fast 
nonadiabatic turning up of the lattice, we study here, 
only about 13-15 wells are occupied after the ramping. 

As the lattice is raised the tunneling amplitude of the 
atoms between the neighboring sites decreases exponen- 
tially and the system becomes more strongly interacting. 
The strong interactions in the lattice enhance quantum 
fluctuations, eventually destroying the long-range phase 
coherence of the atoms. 

The relative phase between the atoms in different lat- 
tice sites during the turning up of the lattice potential 
generally indicates the flow of atoms towards the edges 
of the atom cloud, as the cloud radius increases due to the 
lattice. The dynamics of the phase on a single stochastic 
realization also undergoes stochastic fluctuations due to 
the vacuum noise. The deeper the lattice potential, the 
larger are the variations in the dynamical trajectories 
of the phase between individual stochastic realizations. 
These variations, once averaged over a large number of 
realizations, result in the loss of phase coherence. 

In Fig. 3 we show the normalized phase coherence 
C\ for the atoms between the central well and its first 
neighbor, as defined in Eq. (16) at different initial tem- 
peratures Ti, for two different values of the ramping- up 
time. The phase collapse time exhibits exponential-like 
decrease as a function of 7j for lot = 1. We also show 
the phase coherence between the atoms in the central site 
and in its several nearest neighbors, describing the phase 
coherence along the lattice. The effect of the ramping-up 
time and the final lattice height on the phase coherence 
is shown in Fig. 4. It is interesting to note that in sev- 
eral cases the phase coherence rapidly decays during and 
immediately after the turning up of the lattice poten- 
tial, but remains surprisingly steady at some finite value 
afterwards. 
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FIG. 3: The phase coherence at different temperatures and 
along the optical lattice. We evaluate Ci for the ramping-up 
times lot — 1 (top left; curves from top correspond to the ini- 
tial temperatures kBTi/hw — 3, 9, 15, 21, 27) and lot = 5 (top 
right; curves from top correspond to fcsTi/ftaj = 15,21,27). 
The notch after about one trap period in the coherence func- 
tion is due to a resonance involving the first phonon modes 
Mi and vi and is enhanced with the increasing temperature. 
The estimated phase collapse time r c (bottom left) as a func- 
tion of the initial temperature, obtained from the data in the 
top left diagram by choosing the time when G\ = 0.65. The 
coherence between the atoms in the central site and in the 
first five nearest neighbors (bottom right; curves from top G% 
with i from one to five) for ksTi — 27hio. In all the plots 
N g 1D = nOfkol. 
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FIG. 4: The phase coherence for different ramping-up times 
and lattice height. The same system as in Fig. 3 at ksTi = 
27hu for different ramping-up times (on left; curves from left 
lot = 1,2,3,4,5) and for the final lattice heights s = 37.5 
(lower curve) and s = 25 (upper curve) for lot — 1 (on right). 



C. Optical lattice with a larger number of sites 

In this Section we present the main results of the paper 
by considering the splitting of a harmonically trapped 
BEC by a rapidly varying optical lattice potential. We 
choose the lattice spacing to be d = irl/8 in Eq. (20). We 
then have about 30-35 lattice sites within the classical 
diameter 2R of the BEC. A similar number of sites has 
also been realized in recent experiments in a cigar-shaped 
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FIG. 5: The normalized ground state atom density (thick 
line) in a harmonic (left) and in a combined harmonic and 
optical lattice (right) with s = 20. For illustrative purposes 
we also show the potentials (thin line; in arbitrary units). 
Here Nogm = WOhtol and kl — 8. 



trap with d ~ 2.7/im [4]. The ground state atom number 
Nq = 2000 before turning up the lattice results in the 
occupation number of about hq ~ 90-100 atoms in the 
central site of the optical lattice. The normalized ground 
state atom density is shown in Fig. 5 with (for s = 20) 
and without an optical lattice. 



1. The effect of the final lattice height 

In this Section we vary the final height of the optical 
lattice potential when a harmonically trapped BEC is 
fragmented by means of continuously turning up the lat- 
tice. We evaluate the dynamics of the phase fluctuations 
between the atoms in different lattice sites and the atom 
number fluctuations in individual sites. The system is 
closely related to the recent experiment by Orzel et al. 
[6] where the atom number squeezing in strong optical 
lattices with large filling factors was observed. However, 
the notable difference between our model and the exper- 
imental set-up is that in Ref. [6] the ID optical lattice 
exhibited a very weak radial confinement and did not 
produce a tightly-elongated ID gas with negligible radial 
excitations. 

In Fig. 6 we show the phase coherence between the 
atoms in the central well and in its nearest neighbor C\ 
[Eq. (16)] and the number fluctuations Ano [Eq. (15)] in 
the central well for different final heights of the periodic 
potential at T = 0. The ramping- up time is fixed to lot — 
6 for s = 30, 40 and to lot = 10 for the shallower lattices. 
For shallow lattices the phase coherence remains high and 
steady, but for larger s it is reduced and becomes strongly 
oscillatory. The enhancement of phase fluctuations in 
deeper lattices is associated with progressively increasing 
atom number squeezing. The number squeezing can be 
accurately fitted according to 

^^0.03 + 0.5 e -*/ 8 . (23) 

Due to the large occupation numbers, An are 
strongly sub-Poissonian, approaching the asymptotic 
value (Ano) 2 /no — 0.03 <C 1 for large s. The numer- 
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FIG. 6: The phase coherence between the central well and 
its nearest neighbor Ci as a function of time (top left) at 
T — for different final heights of the optical lattice (curves 
from top represent s = 5, 8, 10, 15, 20, 30, 40 with the atom 
number in the central well no — 90-100). The ramping- up 
time lot = 10, except for s — 30, 40, lot — 6. Top right: Di for 
s = 5, 8, 10 (curves top to bottom). The number fluctuations 
Ano for the central well (bottom right) for the same runs. 
Here g\o = 0.05hiol and Nq = 2000. The number squeezing 
(bottom left) as a function of the final lattice height can be 
accurately fitted according to Eq. (23). 



ics proved more demanding when the final height of the 
optical potential reached 30-40-E r . 

Although the discrete tight-binding Hamiltonian 
Eq. (Al) is only valid for weakly excited and very deep 
lattices, it is interesting to compare the TWA to the 
Bogoliubov approximation of the BETH (Al), introduced 
in Appendix A. The discrete Bogoliubov result for the 
ground state in a uniform lattice [Eq. (A16)] predicts the 
phase fluctuations A<^>oi to be weak when the hopping 
amplitude J times the lattice site occupation number n 
is much larger than the on-site interaction U. From the 
numerical results for D\ [Eq. (17)] in Fig. 6 we can de- 
duce (A(^o,i) 2 and for shallow lattices we find it to be 
roughly by the factor of two smaller than the uniform 
ground state result Eq. (A16). The small phase fluctua- 
tions of the TWA results may be better understood when 
we compare them to the phase fluctuations obtained by 
solving the discrete Bogoliubov theory in a combined har- 
monic trap and an optical lattice; Fig. 14. Including the 
harmonic trap can significantly reduce the phase fluctu- 
ations close to the trap center and enhance them close to 
the edge of the atom cloud. 

The Bogoliubov result of the discrete tight-binding 
Hamiltonian Eq. (Al) for the ground state atom num- 
ber fluctuations [Eq. (A15)] indicates that An* becomes 
strongly sub-Poissonian (An^) 2 -C m, when riiU 3> J. 
The numerical TWA results for the atom number fluctu- 
ations in shallow lattices in Fig. 6 are clearly higher than 
the homogeneous ground state result in Eq. (A15). More- 
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over, in deep lattices (see the discussion in Appendix A), 
the system in the ground state may undergo a quantum 
phase transition from the superfluid to the Mott insula- 
tor state. Here tiiJ <~ U at s ~ 38. However, in the 
simulations we find Auq > 1 for all s. The saturation 
of the atom number squeezing in Eq. (23) and clearly 
higher values of Ano than predicted for the ground state 
indicate that the atoms are not loaded adiabatically into 
the optical lattice. The saturation can be understood if 
we take into account the effects of the finite ramping-up 
time of the lattice. 

In order to preserve the adiabaticity during the 
turning-up of the lattice potential and for the system 
to remain in its ground state, we require that the rate of 
change of any parameter in the Hamiltonian to be slower 
than any characteristic time scale of the system. In opti- 
cal lattices when the lattice height is increased, the most 
rapidly changing parameter typically is the tunneling am- 
plitude for the atoms between neighboring sites. On the 
other hand, at low lattice heights it is more difficult to 
avoid exciting higher vibrational levels within one po- 
tential well, resulting in excitations in the higher energy 
bands. 

The phonon mode energies u>j in the lowest energy 
band decrease with increasing lattice strength [see the 
homogeneous Bogoliubov result for the tight-binding 
Hamiltonian in Eq. (A12)] [60, 61] and, as the lattice be- 
comes deeper, it is progressively more difficult to main- 
tain the adiabaticity with respect to these excitations. 
In deep lattices for the loading to remain adiabatic it is 
therefore required that 



1 dJ(t) 



J{t) dt 



« Wj -(t) 



(24) 



for all the phonon mode energies Uj during the turning 
up of the lattice. 

In Fig. 6 we find the number squeezing to saturate 
around s=20-30, indicating the point where an increas- 
ing number of phonon modes is excited and the phonon 
mode frequencies no longer satisfy the condition Eq. (24). 
Consequently, the turning up of the lattice potential is 
strongly nonadiabatic. Due to the nonadiabaticity, the 
s > 15 cases exhibit significant excess number fluctua- 
tions as compared to the ground state. After a short 
time period over which C\ remains constant, the large 
Arii evolve into large phase fluctuations and C\ becomes 
oscillatory and collapses. 

As argued in Ref . [60] , if the adiabaticity of a phonon 
mode breaks down, the number fluctuations of the mode 
freeze to the value that prevails at the time this occurs, 
i.e., when ujj ~ £(£). If we naively apply this argument to 
the Bogoliubov result in Eq. (A15) for the case where the 
adiabaticity condition Eq. (24) breaks down for several 
modes with uij < ((t), we obtain 



(An) 2 * £ 



2UN V 



where Cj{tj) denotes the value of £(f) when we have u)j ~ 
£(t) for the first time during the turning up of the lattice 
for the jth phonon mode at time tj . Since for all j, Q (tj) 
is in the studied case roughly of the order of u, we have 
the asymptotic value for s — > oo, 



(An,) 2 



Tj 



(26) 



(25) 



qualitatively similar to the results in Fig. 6. 

We may also estimate in this case the time required 
to turn up the optical lattice adiabatically from s = 10 
to s = 35. If, for simplicity, we assume that the optical 
lattice was ramped up in such a way that J(t) ~ Joe~ Kt , 
we obtain in Eq. (24) a constant £ = k. In order to 
have an adiabatic loading of the atoms into the lattice, 
we therefore require that n <C ujj for all the modes j. 
If we assume that the homogeneous Bogoliubov result 
Eq. (A14) for the lowest phonon mode energy oj q , m i n still 
provides at least the correct order of magnitude esti- 
mate at s = 35 with the large filling factor considered 
here, the adiabaticity condition is roughly satisfied if 
k -C 0.1 /u. The corresponding time required to have 
a slow enough ramping from s = 10 to s = 35 is then 
At ^> 40/w ~ 0.6s, where we have used the value of 
the trap frequency uj = 2tt x 10Hz (with the present 
set of parameters this corresponds to the lattice spac- 
ing d ~ 1.3yum for 87 Rb). However, such a slow turning 
up process of possibly tens of seconds in order to reach 
the Mott ground state is extremely difficult to achieve 
experimentally, since the ground state lifetime is limited 
by the losses due to the three-body collisions and the 
spontaneous emission of photons. 

A similar argument also applies to the atom number 
squeezing experiment in Ref. [6]. In Ref. [6] a BEC was 
initially confined in a harmonic trap. The atom cloud was 
then split by continuously turning up an optical lattice 
and the atoms were interfered at different final heights 
of the lattice potential. The atom number squeezing in 
individual lattice sites was calculated from the phase fluc- 
tuations by measuring the loss of visibility in the interfer- 
ence fringes. The lattice had about 12 sites with ^1000 
atoms in the central well. The lattice depths of up to 
s = 50 were reached, corresponding to riiU/J ~ 10 5 . 
Here nJJ was estimated from the size of the atom cloud 
and it varied considerably during the ramping, since the 
lattice light field also changed the transverse confinement 
depending on the lattice height. If we use the same sim- 
plified analysis as before to estimate the time required to 
turn up the lattice adiabatically from s = 10 to s = 50, 
we obtain At » 50ms. Consequently, it is not surprising 
that the ramping-up time 200ms used in the experiments 
was not sufficiently slow to reach the ground state of the 
atoms. The adiabaticity condition At 3> 50ms is not as 
severe as in our simulation example, but in the system 
of Ref. [6] even the lattice height s = 50 is still far below 
the Mott transition point. 

As we already argued in Ref. [30], the requirement of 
a very slow ramping-up time generally makes it difficult 
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to reach the Mott insulator ground state experimentally 
when there are many atoms per lattice site for a large 
number of sites. This indicates that producing techno- 
logically interesting strong atom number squeezing with 
large occupation numbers has serious limitations. In the 
lattices with small filling factors, that have been so far 
used in the Mott insulator state experiments, the above 
argument is no longer valid due to rapidly emerging en- 
ergy gap and much larger quantum fluctuations even in 
a shallow lattice. 

We have shown that for large filling factors the squeez- 
ing of atom number fluctuations within individual lattice 
sites, resulting from the turning up of the lattice poten- 
tial, saturates for deep lattices. The saturation of the 
number squeezing for strong lattices was experimentally 
observed in Ref. [6]. As we already pointed out, such 
a system is not tightly elongated, but we can still make 
qualitative comparisons to the experimental data. Al- 
though the saturation was assumed in Ref. [6] to be an 
artifact of the analysis method of the interference mea- 
surement, we also numerically find similar saturation ef- 
fect in a qualitative agreement with the experiment, with- 
out including the effects of the interference experiment. 
As we argued here, such saturation may result as a con- 
sequence of the nonadiabaticity of the loading process. 

If the loading is sufficiently rapid or the final lattice 
sufficiently high, so that the adiabaticity breaks down for 
a large number of modes, the optimal number squeez- 
ing is proportional to the ramping speed itself and the 
nonlinearity. Both in Fig. 6 and in Ref. [6] the squeez- 
ing saturates at about 15dB when riiU/J ~ 10 4 . The 
ramping- up time r ~ 4Q00H/E r in Ref. [6] is one order 
of magnitude longer than in Fig. 6, but this is compen- 
sated by the weaker hopping amplitude J, so that the 
saturation roughly occurs at the same value of uj n r. 

Due to the nonadiabaticity of the loading of atoms into 
the lattice, the system in Ref. [6] may not have been in 
its ground state during the turning up of the lattice po- 
tential. Similarly to our numerical example, the nonadia- 
baticity in the experiment can induce larger phase fluctu- 
ations than those in the Heisenberg minimum uncertainty 
state for the atom number and phase. Additionally, an 
inhomogeneous phase profile can reduce the visibility of 
interference fringes similarly to phase fluctuations. Con- 
sequently, not all the experimentally detected phase noise 
in the loss of interference fringes may directly correspond 
to the atom number squeezing. Although the increased 
phase noise in the experiment therefore did not provide 
an entirely conclusive measure of the atom number fluc- 
tuations, our numerical simulations, nevertheless, seem 
to indicate that a considerable number squeezing may 
have been present also in Ref. [6]. 

The effect of increased phase fluctuations due to a 
rapid turning up of the lattice may even be more pro- 
nounced in the absence of strong radial confinement, as 
in Ref. [6], as a result of the nonlinear coupling between 
the radial and the axial modes. The numerical solution of 
GPE in ID, in Ref. [31], and in 3D, in Ref. [32], was com- 




FIG. 7: The phase coherences Ci (top left), C5 (bottom left) 
and the number fluctuations Ano (top right) for initial tem- 
peratures (curves from top) fcaTi/fej — 0, 12.5, 22.2, 33.3, 38.5 
(C5 also with 28.5). The phase collapse time t c (bottom right) 
is evaluated at C5 = 0.5, subtracting the ramping-up time. 
Here N g 1D = WOhujl, No = 2000, lit = 10, and s = 20. 

pared to the experimentally observed phase noise using 
the parameters of Ref. [6]. Although this model cannot 
incorporate any quantum effects, the numerics showed 
signs of a reduced interference visibility due to inhomo- 
geneous phase profile [62]. 

2. The effect of initial temperature and nonlinearity 

Finite initial temperature increases the initial noise in 
the Wigner distribution [Eq. (9)] of the TWA as a con- 
sequence of the finite temperature noncondensate frac- 
tion in the excited levels. This is expected to affect the 
coherence between the different sites (as we saw in Sec- 
tion III B) and the fluctuations in the atom number Ano. 
In Fig. 7 we show Ano and the phase coherence G\ for 
different initial temperatures Tj for final lattice height 
a{r) = 20. Here (Ano) 2 increases exponentially as a func- 
tion of Ti. The phase coherence C5 between the central 
well and its fifth neighbor decays significantly faster than 
C\ . The dependence of the phase collapse time r c on the 
initial temperature Tj is approximately linear (compare 
to Fig. 3). At s = 20 the effects of the harmonic trap 
are already significant in C5, since the potential energy 
difference due to the harmonic trap between the central 
lattice site and its fifth nearest neighbor exceeds the tun- 
neling energy V h (j = 5) - V h (j = 0) ~ 2ftw > n Q J. 

In order to estimate the effect of the nonlinearity we 
varied Nogm / (huil) from 100 to 200 for the case of the 
final lattice height s = 20, the ramping-up time lot = 10, 
and the initial temperature = 0. In Fig. 8 we show the 
results for the phase coherence C5 and the atom num- 
ber fluctuations Ano m the central well. The depen- 





FIG. 8: The coherence function C5 (top left) and the number 
fluctuations in the central well (top right). The curves from 
top represent N g 1D /(hujl) = 100, 125, 150, 175, 200. The col- 
lapse time (evaluated when C5 = 0.5 and subtracting the 

ramping-up time) (bottom left) and Ano/n^ 4 (bottom right) 
as a function of the on-site interaction U (in the units of Twj) 
at Tt = 0. Here N = 2000, s = 20, and ujt = 10. 



dence of the phase collapse time on the nonlinearity is 
obtained by choosing the time when C5 = 0.5. We com- 
pare the numerical TWA results for Ario to the analytic 
tight-binding result in the ground state of the uniform 
lattice (A17) by fitting the TWA results according to 
(Ano)/ v /n u " oc U c . Here the exponent c of the on-site 
nonlinearity is the fitting parameter. The TWA simula- 
tions yield c ~ —0.3, as compared to the ground state 
result c = -0.5 of Eq. (A17). In Fig. 8, the (An )/v^0 
plot was obtained by averaging over the time interval 4/w 
after the lattice was turned up. 



3. Validity of the classical GP theory 

The classical mean-field theories, in particular the 
GPE, have been very successful in describing the full 
multi-mode dynamics of weakly-interacting, harmoni- 
cally trapped atomic BECs. However, the GPE has se- 
vere limitations in optical lattices where the interactions 
are enhanced, since the classical GPE disregards thermal 
and quantum fluctuations, decoherence, and the informa- 
tion about quantum statistics. It is especially interest- 
ing to study the limits of validity of the GPE in shallow 
optical lattices. In Ref. [43] it was shown that the ex- 
perimentally observed damping of the dipolar motion of 
atoms in a very shallow lattice with s — 0.25 [11] resulted 
from the large ratio guy/N. If the atom number in the 
simulations was increased, or gm/N decreased, while at 
the same time keeping the chemical potential cx Ngiu 
fixed, the damping rate was reduced exponentially and 
the results approached the classical GPE limit. Here we 
study the turning up of the lattice for much larger filling 
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FIG. 9: The coherence C5 between the central well and the 
fifth neighbor (top left) and the atom number fluctuations in 
the central well Auq (top right) for the case of the final lattice 
height s = 5. The curves from top represent the initial tem- 
peratures k B Ti/hu = 0, 6.67, 12.5, 22.2, 33.3, 38.5. The phase 
coherence C5 for s = 8 (bottom left) and for s = 10 (bottom 
right). The curves from top represent the initial tempera- 
tures k B Ti/fho = 0,22.2,28.6,38.5 and 0,22.2,38.5, respec- 
tively. The nonlinearity NogiD ~ WOhuil, No = 2000, and 
the ramping-up time lot — 10. 



factors and, hence, smaller guy/N. 

In Fig. 9 we show the coherence between the central 
well and its fifth neighbor for different final lattice heights 
and different temperatures and the atom number fluctua- 
tions in the central well for s = 5. The number of atoms 
in the central site no — 90-100. For the case of s = 5 
the phase coherence remains high at low initial temper- 
atures. The number fluctuations at Ti — are weakly 
sub-Poissonian Auq ~ 5.5 < 10 and are increased due 
to the initial finite temperature. However, for s = 8 
and s — 10, the TWA results are notably different from 
the classical GPE dynamics, even at zero temperature. 
For s = 8 and s = 10 the loss of coherence at Ti — 
due to vacuum fluctuations is already comparable to the 
loss of coherence at s = 5 due to thermal fluctuations at 
ksTi = 38.5/10-1. (When the initial thermal population is 
close to ten percent.) For s — 10, the atom number fluc- 
tuations are also more sub-Poissonian Ario ~ 3.9 < 10 at 
Ti = 0; see Fig. 6. 



4- Change of temperature during the splitting 

The optical potential also affects the temperature of 
the atom vapor. If the lattice potential is turned up 
adiabatically, the population of each mode remains con- 
stant and temperature T can change dramatically, as the 
contribution of each mode to T changes by the ratio of 
the final and initial mode energies or. /cjj* . An adia- 
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5iE) tlal — to a final value <7^ al during the time r g , be- 
fore starting to ramp up the optical potential. Exper- 
imentally, such a technique could possibly be employed 
by means of using Feshbach resonances to tune the value 
of the scattering length a. 

The field operator for the initial state of the noninter- 
acting system is that of the ideal harmonic oscillators in 
thermal equilibrium: 



FIG. 10: The contribution of the lowest modes to temperature 
(left) for ksTi/hio — 12.5 (dashed line) and 37 (solid line). 
Curves from top to bottom represent ujt — 10, 20, 30 in both 
cases. The average temperature of the first five modes after 
the ramping ujt = 30 (right) for an initial temperature fc_sTi = 
0, 12.5,22.2,37. Here 3115 = 0.015/kjZ, iV = 2000, and s = 5. 



batic increase in the lattice strength may both increase 
or lower T, depending on whether the excited band is 
occupied [63] and in the experiments the condensation 
temperature has been found to be sensitive to the lattice 
height [64]. In Fig. 10 we estimated the 'temperature' 
of several lowest phonon modes in the TWA simulations 
by evaluating the corresponding occupation numbers rik . 
This was obtained by calculating the projection of ipw 
to the Bogoliubov modes of the BHH using Eq. (A6), as 
explained in Appendix A. The averages are taken over 
a time period before any significant rethermalization oc- 
curs after the ramping [65] . The modes 2 and 4 are highly 
excited for the case of short r, due to the nonadiabatic 
loading. The excitations are damped out at higher initial 
temperatures Tj and for the case of the slowest ramping 
ujt = 30, representing the situation where uj-it 1 and 
0J4.T 1. It is interesting to note that the excitations 
of the forth mode are only damped out when the rate of 
change in the tunneling amplitude £ is much smaller than 
the corresponding mode energy, or when uj<± ~ 26£(r). 
This is more restrictive condition than the one found in 
Ref. [60]. For very fast ramping ujt < 3, the variation 
of Ti is already completely dominated by the excitations 
due to the rapid turning up of the lattice. 



^* = °)=^T7I ex P ("^) 

+ ajRjHjix/l) exp (-^) , (27) 

where we have explicitly separated the BEC mode, Hj 
is the jth Hermite polynomial, and Rj = (y^ 3 j!/) -1 / 2 . 
The sampling of the Wigner distribution for aj and ctQ 
to generate ipw(x, t = 0) is then performed exactly as in 
the interacting case. 

In Fig. 11 we show the TWA simulation results for the 
initially noninteracting system for which we first slowly 
turn up the interactions to the value gm — 0-05HujI, be- 
fore ramping up the optical lattice potential. We also 
show the results for the initially interacting system for 
which the initial state is obtained by solving the Bogoli- 
ubov equations with g\r> = O.Obhujl. The initial temper- 
ature of the noninteracting case is set to be equal to the 
temperature of the interacting system we want to com- 
pare. Therefore the condensate and the noncondensate 
populations in the initial state of the TWA simulations 
are slightly different in the interacting and noninteracting 
cases. 

In the example case the two approaches yield surpris- 
ingly similar results at low Ti. At higher Ti the atom 
number fluctuations in the case of the ideal gas initial 
state are smaller. This is most likely due to the lack of 
any notable rethermalization in the ID TWA simulations 
[65] , so the turning up of the nonlinearity fails to produce 
the thermal equilibrium state. 



5. Turning up the nonlinearity 

The numerical solutions for the initial equilibrium 
state in the TWA simulations may in some cases, es- 
pecially in higher dimensions, be difficult to obtain. In 
Ref. [42] the ideal, noninteracting condensate was used as 
an initial state for the TWA simulations, but the inter- 
action constant was first continuously ramped from zero 
up to some final value, before the actual atom dynamics 
was studied. By means of turning the nonlinearity up 
slowly enough, the goal was to produce the initial state 
of an interacting system. If the ramping is slow enough, 
the dynamics is then expected to resemble to that of the 
interacting system. Here we study the atom dynamics 
by first linearly turning up the interaction constant from 



6. Sampling of the noise for the initial ground state mode 

In the TWA simulations the quantum fluctuations of 
the initial BEC mode are included by assuming the BEC 
to be in a coherent state, so that the initial state is sam- 
pled according to the corresponding Wigner distribution 
in Eq. (11). Although the condensate mode is expected 
to exhibit sub-Poisson atom number fluctuations, the ad- 
vantage of the coherent state representation is that the 
Wigner distribution is positive. For instance, the Wigner 
distribution of a number state with n atoms is nonpos- 
itive: W(a,a*) = 2(-l)" exp (-2|a| 2 )L„(4|a| 2 )/7r [48], 
where L n denotes the Laguerre polynomial, and cannot 
be represented in the TWA simulations without doubling 
the dimension of the phase space [66] analogously to the 
positive P representation [47]. However, as we already 
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FIG. 11: The phase coherence and the number fluctuations for 
an initially noninteracting system (c/i£> tlal = 0) with the inter- 
actions linearly turned up to the value gilf 1 = giD = 0.05hujl 
(dashed red line) and in an initially interacting system (solid 
black line), obtained by solving the Bogoliubov equations for 
gio = 0.05hujl with iV = 2000, s = 20, and ujt = 10. The 
initial temperature is knTi/hw — 6.67 (on top) and 33.33 
(at bottom) and the ramping-up time of the nonlinearity 
uiTg — 15. The results for the evolution of the ideal gas 
describe its dynamics after the end of the ramping of the 
nonlinearity. 



noted earlier, the particular representation is not so cru- 
cial, since the main contribution to the matter wave co- 
herence is due to the thermal and quantum fluctuations 
of low-energy phonons and the quantum fluctuations of 
the initial state of the BEC mode are not very important. 
In Fig. 12 we show the TWA results for the typical physi- 
cal parameters of our system both when the BEC mode is 
sampled according to Eq. (11) (with Poisson atom num- 
ber fluctuations) and when we entirely ignore the quan- 
tum fluctuations of the BEC mode and treat it classically 
having a fixed number of atoms. A classical treatment 
of the ground state does not affect the prediction for the 
phase coherence along the lattice, but produces slightly 
smaller atom number fluctuations in the lattice sites. 



7. Importance of the symmetric ordering 

The TWA simulations return quantum expectations 
values for the operators that are symmetrically ordered. 
This generally creates significant complications in inter- 
preting the simulation results in terms of normally or- 
dered expectation values, as we already emphasized in 
Section II B. However, as we here demonstrate, the cor- 
rect transformation between the different operator or- 
derings is very crucial in obtaining the correct physical 
results. 

Since the P-representation of the density matrix re- 
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FIG. 12: The phase coherence Ci and the atom number fluc- 
tuations Ano for cases where the initial ground state was sim- 
ulated classically with a fixed number of atoms (solid black 
line) and where the initial ground state included quantum 
fluctuations while assuming a coherent state (dashed red line). 
Here % = 0, N g 1D = lOOfojZ, iV = 2000, s = 20, and 

UJT = 10. 

turns the expectation values of normally ordered oper- 
ators, e.g., (a*ctj)p — (a^dij), we may also call the 
stochastic dynamics that assumes the operator expecta- 
tion values to be normally ordered as the "truncated P 
approximation" (TPA) [38] . The P-representation is sin- 
gular at T = 0, but is often used to study finite temper- 
ature systems. A review of its mathematical properties 
may be found in Refs. [47, 48]. 

We illustrate with a numerical example the importance 
of the particular stochastic representation of the density 
matrix. We integrate both TWA and TPA according to 
Eq. (5) for the same system at finite initial tempera- 
ture Ti by means sampling the Wigner and the P dis- 
tributions. The results shown in Fig. 13 for the phase 
coherence C\ and the atom number fluctuations Ano are 
then useful in understanding the importance of the oper- 
ator ordering. The P-distribution does not introduce the 
correct noise in the initial state and therefore it fails to 
produce the right phase coherence properties and under- 
estimates the atom number fluctuations. At high tem- 
perature, the TPA is closer to the TWA because thermal 
noise starts dominating vacuum fluctuations. 

IV. CONCLUDING REMARKS 

We studied the loading of a harmonically trapped BEC 
into an optical lattice. As the lattice height is increased, 
quantum fluctuations become more important and the 
classical GP description breaks down. We found that 
the TWA provides a powerful tool to study nonequilib- 
rium dynamics of bosonic atoms in an optical lattice. We 
showed that, especially at low temperatures, the correct 
quantum statistics of phonon modes and the accurate 
treatment of the operator orderings in the Wigner repre- 
sentation are very important. 

The TWA incorporates the full multi-band description 
of the lattice dynamics and becomes more accurate when 
the occupation number of the sites is large and quan- 
tum fluctuations are not too dominating. In the TWA 
simulations of the turning up of the lattice potential we 
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FIG. 13: The phase coherence and the atom number fluctua- 
tions, obtained by taking into account the symmetric opera- 
tor ordering in the TWA (dashed line) and by assuming the 
stochastic representation to be normally ordered (solid line) 
for an initial temperature of the system ksTi/hu! = 6.67 (left) 
and 38.46 (right). The same system as in Fig. 12. 



APPENDIX A: BOGOLIUBOV THEORY IN A 
DISCRETE TIGHT-BINDING APPROXIMATION 

It is useful to compare the TWA simulation results 
to the ground state fluctuations of the atoms, obtained 
from the Bogoliubov theory. In a deep optical lattice, 
with very large s, and close to the ground state only one 
mode per lattice site is important and the system can 
be approximated by the discrete Bose-Hubbard Hamilto- 
nian (BHH), where we expand the field operator on the 
basis of the Wannier functions and keep only the lowest 
vibrational states in each lattice site, ip{x) — J^- bif]i{x): 

H = J2 fc&ft " J(btk+i + H.c.) + ^(b\) 2 b 2 ] , (Al) 



where the summation is over the lattice sites, J ~ 
— J dxr/* (x)Ci]i + i (x) is the hopping amplitude between 
the nearest-neighbor sites, U ~ gm J dx\r)j(x)\ is the 
on-site interaction constant, and v.j = j 2 d 2 muj 2 /2 repre- 
sents the harmonic trapping potential, with j — site at 
the trap center. We may approximate the Wannier func- 
tions r\i by the ground state harmonic oscillator wave 
function with the frequency 



2s 1 ' 2 E r 



(A2) 



found the atom number squeezing to saturate for deep 
lattices, which can be explained by the finite turning-up 
time. It is numerically more demanding to study sig- 
nificantly longer ramping-up times that would result in 
considerably more reduced atom number fluctuations. It 
is also expected that the TWA would eventually break 
down over longer time scales very close to the MI state 
that could be interesting to investigate. In our theo- 
retical study of dissipative dipole oscillations of bosonic 
atoms [43] , we successfully modeled the system with the 
TWA for the case of considerably stronger quantum fluc- 
tuations than in the present work, in a good agreement 
with the experimental results [11]. This seems to suggest 
that also notably stronger atom number squeezing could 
be studied within the TWA than the one reported in this 
paper by considering very long ramping-up times. As our 
analysis shows, however, that, in the case of lattices with 
large filling factors, very long ramping-up times, that are 
required to reach the MI ground state, can be extremely 
challenging in actual experiments. 



The TWA simulations could also be extended to 2D 
or 3D lattices. In 2D and 3D the implementation of the 
TWA may require a special care, but, at least in the uni- 
form space, it has been successful in modeling thermal 
fluctuations [38]. In higher dimensions also more com- 
plex initial states could be considered, such as topologi- 
cal defects and textures that may be prepared by phase 
imprinting [67]. 



which is obtained by expanding the optical potential at 
the lattice site minimum [68]. Analytic approximations 
for U and J may be derived using the Gaussian approx- 
imations to the Wannier wave functions or the Matthieu 
functions of the energy band. When we compare the 
TWA results to the BHH, we frequently extract the ex- 
pectation values involving b using Eq. (13) with b ~ a. 
We tested that using different projections does not affect 
the results. For mJ > U, with n,; = (b\bi), the system 
is in the superfluid regime with the long-range phase co- 
herence and is expected to undergo the MI transition 
at 2.2riiJ ~ U [69, 70], resulting in a highly number 
squeezed ground state. 

In the Bogoliubov approach to the BHH we introduce 
the discrete position coordinate is x — jd, where j is 
the integer labeling the lattice sites. We decompose the 
ground state and the excited state fractions into mutu- 
ally orthogonal subspaces and study the linearized fluc- 



tuations around c,- 



h = c j 



Sbj 

Yi\fMxn-Kijd}xi], (A3) 



where Xn are the quasiparticle operators and [/ n (j<2), 
h n (jd)] are the Bogoliubov modes, obtained from: 

£fti + Aft 1 + ft 1 ) - <? s uk = e 3 p n 

L d b? n + J(h{ +1 + h^ 1 ) - cfU.fi = -esK , (A4) 
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where Cj = 2njU + i/j — u and /■£ = f n (jd). 

The expression for the fluctuations 5b j in Eq. (A3) may 
be inverted by using the relation: 



suitable for an analytic treatment. In the case of the ho- 
mogeneous system with n atoms per site, we obtain the 
excitation frequencies [60, 61]: 



E [U3d)r m {3d) - h n (jd)h* m (jd)} = S mn , (A5) (^) 2 ^ 4Jsin2 



4 J sin 



2nU 



(A12) 



in order to obtain x„: 

= E [fnU^Sbj + h n (jd)5b] 



(A6) 



In Section III C 4 we use this expression to evaluate the 
occupation numbers 



n k = (XfeXfe) 



(A7) 



of the low energy quasiparticlc modes in the lattice in the 
TWA simulations by substituting 5b j = bj — Cj, so that 
Hk is obtained from the numerically evaluated correlation 
functions involving b/s. 

The atom number operator hi at the site i may be 
obtained by noting that we may expand to first order in 
5b f- b]h ~ rij + n,, for 



(A8) 



where Wj = fj — hj. Moreover, we introduce the phase 
operator at the site i as 



<Pi 



2,/n, 



(A9) 



for which the commutator [hj, <j>j] = i and we have de- 
fined Tj = fj + hj . 

Diagonalizing Eqs. (A4) yields the results for the num- 
ber fluctuations in the zth site and the phase fluctuations 
between the kth and Zth site: 



(An,) 2 = n l J2\w 3 \ 2 {2h :j + 1) . 



(A10) 



and 

(Atpki) 2 



{(<Pk - <f>if) 



Tj{kd) Tj{ld) 



(2n j + l), (All) 



where hj = (XjXj) is the thermal population of the jth 
phonon mode in the lattice. 

In an inhomogeneous BEC, the phonon modes are spa- 
tially localized and Eqs. (A10) and (All) generally lead 
to spatially varying density and phase fluctuations. How- 
ever, in the homogeneous case (yi — 0) the Bogoliubov 
expressions for the phase and number fluctuations are 



where q is the quasiparticle momentum with periodic 
boundary conditions, so that qd = 2-Km/N Pl for m — 
— N p /2, . . . , N p /2 — 1, and N p denotes the number of lat- 
tice sites. Moreover, in the experimentally interesting 
regime nlJ > Jwe obtain 



hhjq ~ 2V2nJU |sin (qd/2)\ , 
and for N p ^> 1 the lowest phonon mode energy 

2irV2n~JU 



huj, 



q.mm 



N r , 



(A13) 



(A14) 



Moreover, in the limit nil » J, the expressions for the 
number and phase fluctuations in Eqs. (A10) and (All) 
can be approximated by [60, 61]: 



(An,) 2 -E 



4nJN r 



■(2n, + l), 



(A15) 
(A16) 



where h q is the occupation number of the phonon mode 
with the quasimomentum q. By replacing the sum by an 
integral when N p is large, these yield at T = 0: 



(An, 



1 /8nJ 
"tt V~ 
1 (2U 



1/2 



2 V nJ 



1/2 



(A17) 
(A18) 



Formulas (A17) and (A18) provide useful comparisons to 
our numerical results. 

In Fig. 14 we show the numerical results of the number 
and phase fluctuations in a harmonic trap, obtained by 
solving the Bogoliubov approximation to the BHH from 
Eqs. (A4) for a typical system in our simulations. We 
evaluate Ario using Eq. (A10) and compare it to the ho- 
mogeneous results with n = no, determined by Eq. (A17). 
We also calculate (At^i) 2 as a function of the lattice 
height s using Eq. (All) for several sites along the lat- 
tice and we evaluate (A^?o,i) 2 in the homogeneous case 
fromEq. (A18). The inhomogeneous results for the phase 
fluctuations are numerically very unstable due to the 
summation formula in Eq. (All) and spatially localized 
phonon modes. Consequently, the results for (A<p ,i) 2 
should only be considered as an order of magnitude esti- 
mates, describing the qualitative behavior. Numerically 
we find the Bogoliubov result in a harmonic trap for Ario 
to be larger and for A<p i smaller than the homogeneous 
result. The harmonic trap significantly reduces phase 
fluctuations close to the trap center and enhance them 
close to the edge of the atom cloud. 
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FIG. 14: The number and phase fluctuations of atoms in a 
combined harmonic trap and optical lattice, obtained by nu- 
merically solving the Bogoliubov modes in the discrete tight- 
binding approximation. On left we show Ano for various val- 
ues of s (solid line) and the analytic homogenous result with 
n — no in Eq. (A17) (dashed line). Here the total atom num- 
ber No = 2000, the nonlinear constant g\r> = 0.05fkjl, T — 0, 
and no —80-90 atoms. On right for the same inhomogeneous 
system we show approximate (Atpo t i) 2 for i = 1, 13, 16 (solid 
lines from bottom) and the analytic homogeneous result for 
(A<£>o,i) 2 (dashed line) from Eq. (A18). The inhomogeneous 
results for (Aipo,i) 2 should only be understood qualitatively 
as an order of magnitude estimates. 



APPENDIX B: THREE-BODY LOSSES 

The TWA approach in this paper ignored any atom 
losses due to three-body collisions. The TWA with in- 
corporated three-body losses would result in a dynami- 
cal stochastic noise term in Eq. (5) and more demanding 
numerics. We may estimate the importance of the three- 
body losses in an optical lattice system as in Ref. [71]. 
In a deep lattice with a fragmented condensate the atom 



loss rate at the central lattice site (with tiq ^> 1) is ap- 
proximately 



dno 



-Tnl 



12 



d 3 r\r, Q (r)\ ( - 



(Bl) 



where K% is the recombination event rate and 770 ( r ) is 
the ground state wave function in the central site with 
the atom number no- For simplicity, we have ignored the 
correlations between different lattice sites. 
We obtain from Eq. (Bl): 



A', 



(9idY 



67A- 3 
a 2 / 4 



(B2) 



Here on the right-hand side we have introduced the same 
parameters as in Fig. 6 with s = 20. In order to the 
three-body loss rate to be negligible over the time scale 
of the simulations we require Tn^ <SC lo that may be 
satisfied for sufficiently large I. For the ID description 
of the dynamics to be valid we also need to maintain 
huj± > nU ~ ngmks 1 / 4 /v2~7r. For example, for 87 Rb we 
use K 3 ~ 2.2 x 1CT 28 cm 6 /s [72] and a ~ 5.313nm. Then 
u> = 2tt x 1Hz yields Tuq ~ 0.006w and, using an aver- 
age occupation number n, hu± ~ Zngrnks 1 ' ' 4 / 'v2tt with 
lu± 3> uj. The effect of the three-body losses could be fur- 
ther reduced by using different atoms and the Feshbach 
resonances. 
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